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Abstract 

Boundary integral methods for the solution of boundary value PDEs are 
an alternative to 'interior' methods, such as finite difference and finite element 
methods. They are attractive on domains with corners, particularly when the 
solution has singularities at these corners. In these cases, interior methods can 
become excessively expensive, as they require a finely discretised 2D mesh in 
the vicinity of corners, whilst boundary integral methods typically require a 
mesh discretised in only one dimension, that of arc length. 

Consider the Dirichlet problem. Traditional boundary integral methods ap- 
plied to problems with corner singularities involve a (real) boundary integral 
equation with a kernel containing a logarithmic singularity. This is both te- 
dious to code and computationally inefficient. The CBIEM is different in that 
it involves a complex boundary integral equation with a smooth kernel. The 
boundary integral equation is approximated using a collocation technique, and 
the interior solution is then approximated using a discretisation of Cauchy's 
integral formula, combined with singularity subtraction. 

A high order quadrature rule is required for the solution of the integral 
equation. Typical corner singularities are of square root type, and a 'geometri- 
cally graded h-p' composite quadrature rule is used. This yields efficient, high 
order solution of the integral equation, and thence the Dirichlet problem. 

Implementation and experimental results in matlab code are presented. 
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1 Introduction 



This report describes a research project carried out from March to October 1992, at 
the Department of Mathematics, The University of Queensland, Australia. It was 
carried out under the supervision of Dr Graeme A. Chandler, and was accredited 
as a #30 project, coded MN882. 

Techniques related to the CBIEM have been analysed in , and used in ||] . 
The CBIEM is also closely related to the 'Complex Variable Boundary Element 
Method' p5| . This report contains an application of it, using h-p quadrature to 
achieve high rates of convergence, even in the presence of corner singularities. This 
application owes its conception to my supervisor. 

The ideas of graded meshes and h-p quadrature (numerical integration) are pre- 
sented in §eI and are illustrated by experimental results. The CBIEM itself is de- 
scribed in §p|. §^ details numerical implementation of the CBIEM, using the quadra- 
ture technique described in and presents error results for some test problems. 

concludes the report with suggestions for further development, matlab code 
written for the implementation is listed in Appendix 
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2 h-p Quadrature Methods 



2.1 Introduction 

This section describes a high order numerical integration (quadrature) technique, 
that retains its high order in the case of end point singularities in the integrand. The 
method uses a graded mesh, with integration rules of high order used on larger inter- 
vals, and low order on smaller intervals. To achieve the 'best' possible convergence 
rates, whilst including the end points of each interval, the basic quadrature rules 
used are Gaufi-Lobatto. The underlying mesh is graded in a geometric manner. As 
the method of using different quadrature rules on internal intervals is a generalisa- 
tion of earlier 'h' and 'p' methods, the resulting composite quadrature rule is called 
a 'geometrically graded h-p^ method Q. 



2.2 Quadrature Methods — the Questions 

Given an integrand / : [a, b] M, consider the numerical approximation of the 
definite integral by a rule {x^, Wk} on n points: 

rb 



(Xk) Wk 

k=l 



The interval [a, b] is possibly infinite or semi-infinite, but this report considers only 
finite intervals; and without loss of generality, let [a,b] — [0, f]. Similarly, the inte- 
grand could include a weighting factor uj (x), but this is not required here. 

The degree of a quadrature rule is the maximal degree of the polynomial that it 
can integrate exactly]^ That is, if the degree of a rule on n points is p, then: 



x^dx = ^ xlwk, j = -.p. 

k=l 



If / is smooth, the rate of convergence for n point Gaufiian quadrature is O (p") 
(for some p < 1), and for the composite Simpson's rule it is O {n^'^)- That is, the 
error decreases more quickly for GauBian quadrature. This is not true in general 
if / has a singularity.^ For example, consider the 'square root' singularity f{x) = 
\/x, X e [0, 1]. In this case, the rate of convergence for GauBian quadrature falls 
to O (n^'^) , whilst that of Simpson's rule is O (n^'^/^) . Even so, using a composite 
Simpson's rule and a graded mesh, a convergence rate of O {n~'^) can be recovered. 

A composite quadrature rule is created by subdividing the interval of integration 
into m subintervals, and evaluating the integral over each subinterval using an 
appropriate quadrature rule. Choosing Xj-i < Xj, j ~ 1 : m, and xq — 0, Xm = 1: 

/ f{x)dx^Y. / fi^)d^- 
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^Comments on errors refer to discretisation, not machine roundoff error unless explicitly stated. 

■^'Singularity' is intended to always mean 'end point singularity'. If a particular singularity in 
the integrand is not at an end point, then the interval can be subdivided so that the singularity is 
at the end points of the two subintervals. Most quadrature methods perform poorly on integrands 
with internal singularities. 
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Grading the mesh means that the subdivision is organised in some way. A de- 
scription of a generahsed mesh grading to cater for compUcated possibihties, such as 
adaptive quadrature is found in Q. Here, the simphfication of nonadaptive meshes 
is used. Meshes are graded by assigning mesh points according to some simple func- 
tion. A quadrature rule of degree pj (a function of rij , the number of points used 
by the rule) is used on the interval [xj^i,Xj]. This raises two issues: 

• How should Xj and Uj be chosen? That is, how should the mesh be graded, 
and how should the degree of the quadrature rules on each subinterval vary? 

• How is this procedure dependent on the integrand? Consider the generalisation 
of ^/x to \x\" , or more pathological cases. For experimental work, a known 
integrand is easy to deal with. What of more general cases, where no explicit 
functional information is known? 

The remainder of this section describes some partial answers to these questions, and 
displays some experiments. The answers deal with the case — 1 < a < 1, and 
the experiments demonstrate the case a = 1/2. 



2.3 Gaufiian Quadrature 

Gaul3ian quadrature rules are the best possible rules in the sense that they are of 
maximal degree.^ This is due to the exploitation of the maximal number of degrees 
of freedom in the choice of their nodes and weights. GauB methods divide into 
categories depending on the associated weight function, and whether there are any 
prescribed quadrature points. For the applications in this report, it is preferable to 
use the end points of the interval as quadrature points, and a unit weight function 
is assumed. The appropriate set of rules are called the Gaufi-Lobatto rules pages 
101-104]. 

Theorem 2.1 (Gaufi Lobatto Quadrature) 

Given f g (72T1-2 ^ point Gaufi-Lobatto quadrature rule (n ^ 2), has nodes 

a = xi < X2 < ■ ■ ■ < Xn-i < Xn = h, and positive weights wi, . . . ,Wn such that: 

f {x)dx = '^f (Xk) Wk + En. 

fc=i 

Here is dependent on f , a, b and n: 

(2n-l)[(2n-2)!]^^ ^ ' ^ ^ ' ^ ^ ^ 

The rule is of degree p — 2n^S (this is always odd). Observe that forn = 2 the rule 
is the trapezoidal rule, and for n = 3 it is Simpson's rule. Only for n ^ A do these 
rules diverge from the series of closed Newton-Cotes rules (see Table |^ on page |^. 

An algorithm [pr] for finding {xk^Wk} using a matrix eigenvalue technique is im- 



plemented in Appendix A.£ 



^This may of course, not be ideal for the particular application, but in the absence of theoretical 
functional information about the integrand, nothing beyond this can be said about the convergence 
rates of any quadrature method. In practice, with commonly occurring functions, there is a certain 
amount of implicit theoretical information which can be used in error analysis. 
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2.4 Graded Meshes 



A number of different methods for grading meshes appear in the hterature. Three 
important methods |^ are described below. In each case, the mesh subdivides 
the interval [0, 1], when the integrand has a singularity at 0. The meshes have m 
subintervals, that is to + 1 points, including the ends. The jth mesh point is at Xj, 
and the jth interval is of width hj: 

1. Quasiuniform. The mesh is essentially uniform; that is for some constant r < 
1, hj S [ft.T, h] , j = 1 ; TO — 1, where h — maxj {hj}. 



3. Geometric. For some < a < 1, Xj — cr™ j = 1 : to; and xq — 0. 
A geometrically graded mesh is illustrated on one segment of a closed contour 



2.5 h, p and h-p Quadrature Methods 

The '/i-p' nomenclature presented here originated in papers by Babuska et al. 
p^ , on finite element methods, based on previous work which did not explicitly 
use this schema. The following discussion of the three methods refers to their use 
with graded meshes. 

2.5.1 h Methods 

An h quadrature method is composed using two steps: 

1. Choose an underlying mesh of subintervals; possibly a graded mesh determined 
by the user, from analysis of the singularities of the integrand. 

2. Integrate over each of the to mesh intervals, applying the same n point quadra- 
ture rule. The result is a composite quadrature rule on a total of N points. 
These TV points will be called the node points from now on. The functional 
relationship N (to, n) is dependent on whether the basic quadrature rule is 
open or closed. (If the rule is open, the original mesh points are not included 
in the final rule.) Observe that the user cannot arbitrarily select N, only to 
and n. 



A particular basic rule is decided upon (e.g. Simpson's rule), and desired accuracy 
is hopefully attained by simply increasing to (that is, N). Whatever grading is 
chosen, the separation of the node points (h) decreases as to is increased, hence the 
name '/i method'. For a uniform mesh (which works well for smooth integrands), 
h = {b — a) / {N — 1) is constant. Alternatively, open rules, or Gauf5 rules can be 
used, the only important factor is that all the basic rules are of the same type and 
degree. 



2 



Algebraic. For some j ^ I, Xj 




in Figure ^ on page 19 




TO (ri — 2) + TO + 1 (closed basic rule) 
mn (open basic rule). 
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2.5.2 p Methods 

In a p method, again a graded mesh is created. Integration is performed over each 
mesh interval using a basic quadrature rule on n points. Here, n instead of m is 
varied by the user. That is, for a given number of mesh subintervals, m, a set of 
rules of increasing degree (that is n, the number of points involved) is used, until 
desired accuracy is obtained. The same functional relationship N (m, n) exists. As 
n is increased, the rules used grow in their degree (p), hence the name 'p method' 
(see also Table |^) . 

To illustrate, consider the family of closed Newton-Cotes rules. Assume that 
the interval has been subdivided, possibly using an adaptive algorithm that chooses 
smaller subdivisions where there the integrand has greater derivative. Approximate 
the integral over each division using the trapezoidal rule {p — 1), and inspect the 
result. If it is unacceptable, repeat using Simpson's rule {p — 3). Continue this 
process until results are acceptable. 

2.5.3 h-p Methods 

The h-p method is the natural combination of the two previous methods. The 
user may vary both m and n. The idea behind this is to create a composite rule 
that minimises errors in the approximation, for a given number of node points 
N. (Experiment demonstrates that this is achievable.) The user chooses a family 
of basic quadrature rules, then decides how to vary n with mesh interval. As the 
singularities considered are always at end points, a good choice is to organise small 
mesh intervals and low degree rules (small n) near the end points, and larger mesh 
intervals and high degree rules away from them, where the integrand is expected to 
be smooth. 

A simple choice is to begin with a rule on n = 2 points on the smallest interval, 
then linearly increase n with the number of the mesh interval. Other discrete integer 
functions rij, j = 1 : m are easily designed. The only constraint on these functions 
is that if any error analysis is to be done, there should be some regularity in Uj. 
(Choosing basic rules from the same family facilitates this.) This implementation 
uses the Gaufi-Lobatto rule of degree 1 (n = 2) on the first interval, degree 3 (n = 3) 
on the second, etc. 

Creating the composite quadrature rule is quite difficult. Each of the basic 
quadrature rules must be appropriately scaled and shifted, and then coincident 
mesh points must be combined. This is further complicated in the cases of closed 
meshes, closed quadrature rules, and contour integration, where the end points of 
various segments of the parameterisation must also be combined. (This is exacer- 
bated if the contour is closed.) The CBIEM requires all of these to be implemented. 
The (closed) contours involved have corners, and the integrand will usually have 
singularities at these corners. As it will be important to keep the collocation points 
(see § ^.2.3| ) between, and not on, the corners, the underlying meshes must include 
end points. This means that the basic quadrature rules must be closed, so as to 
include the end points. 

The literature recommends using a geometrically graded mesh, with an h-p 
quadrature method. (This is implemented in the CBIEM.) For maximum efficacy, 
the basic quadrature rules chosen must be of maximal degree, which restricts them 
to Gaufi rules. They must also be closed. An n point rule already has two of its 
points fixed, at the ends. The appropriate rule is known as the GauB-Lobatto rule, 
which is of degree p — 2n — 3. 
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2.6 Error Analysis for the h and h-p Methods 



This section is tedious, and consists mainly of technical arguments. The important 



parts are Theorem 2.2 on page K Theorem p.3| on page 03; and the experimental 



results in §§2.6.3 and 2.6.4. The rest can be skipped without loss of continuity. 



2.6.1 Error Analysis for the h Method 

This section computes an error bound for the h method using an algebraic mesh, 
for the integrand |a;|", on [0, 1]. Clearly a > —1 is necessary for the integral to be 
proper, and thus make its computation sensible. For — 1 < a < 0, ^ C" [0, 1], 
so the integrand is unbounded, but the integral is nonetheless defined. If ^ a < 1, 

e C° [0, 1], but |a;|" ^ [0, 1]. If a € N, then £ [0, 1], so the singularity 
vanishes and the case is of lesser interest. If a > 1, but a ^ N, then all of the higher 
derivatives at x = will not exist. Using the notation [xj and [a;] for the least 
integers (respectively) greater than and less than x € K, in general, for a > 1, 

€ CL"J [0, 1], but |a;|" ^ C^"! [0, 1]. These cases are not particularly interesting, 
so the limit a < 1 is made for simplicity. That is, consider — 1 < a < 1, which 
includes the paradigm example x^/^. The nth derivative of / (x) = \x\" , for x 7^ 0, 



(„) _ r{a + i) 



r> (^) = 



r (a + 1 - 71) ' 



Consider the interval [0, 1] partitioned into m subintervals, where Xj is the jth 
mesh point, j = : m, and hj the width of the jth interval, for 

j = 1 : m. Recall, for an algebraic grading, a real constant 7 ^ 1 is chosen,^ and 
the mesh is defined by: xj = {j /m)p , j = : m. Differentiating and applying the 
mean value theorem shows: 

\m J \ m J m< m\mj d] 

For a geometric grading, for some constant < a < 1^ Xj — ct™^^, j — 1 : m; 
xq = 0, so: 

hj = - = (1 - a) (t™-^ j = 2 : m, hi = a™-^ 

Consider an algebraic grading, using closed basic quadrature rules. For an h or h-p 
method, the integral on the jth interval is computed using a quadrature rule with 
rij ^ 2 points. For an h method, Uj is constant; for instance, rij = 2 means that the 
integral over each mesh interval is computed using the trapezoidal rule - the rule 
is a composite trapezoidal rule. The degrees of some common quadrature rules are 
presented in Table 0. 

Now consider the global error of the difference between the true solution / 
and the approximation /,„, induced by the quadrature on m subintervals, that is 
I = Im + Em- For the function |x|", bounds on are readily found. Define Cj 
as the component of due to the jth mesh interval, that is Em — SjLi ^j- To 
bound Em, note Em ^ X^i'li then bound each of the \ej\. 



^Choosing 7 < 1 results in some mesh points possibly lying outside the interval [0, 1]. 
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Degree 


Gaufi-Lobatto 


Newton- Cotes 


1 


2 


2 (trapezoidal) 


3 


3 


3 (Simpson's), 4 (3/8ths) 


5 


4 


5 (Boole's), 6 


7 


5 


7, 8 



Table 1: The number of points n, and associated degrees of some closed quadrature 
rules. For GauB-Lobatto rules, the degree is 2n — 3, and for Newton-Cotes rules, 
the degree is n if n is odd, else it is n — 1. 



The error result for a degree p quadrature rule on an interval [ ] , of width 

hj, with a function / (x) G C^^^ [^j-i: is: 



C(p)/ifV<f+i)(C), 



(2) 



For n point GauB-Lobatto quadrature, the degree is p = 2n — 3. (^ is derived from 
(0) by a scaling argument, and: 



C{P) = - 



(p + 3)(p + 5) [((p-l)/2)!] 
2^p + 2) [{p+l)\f 



(3) 



As g C°° (0, 1], (H) holds for every mesh interval except the first, where the 
error is known exactly, e.g. for the trapezoidal rule: 



1 1 



a + 1 2 



(4) 



As the maximum value of the (p + l)th derivative of |a;|" on the jth interval, is at 
its left hand end, {£,) < (x^-i), the total error can be bounded: 



Here, C refers to a positive constant, independent of to, that may vary from line 
to line. Substituting for hj+i and /(p+i) [xj) using an algebraically graded mesh 
yields: 



m — 1 

E-m < |ei|+C^ 

m — 1 

^ |ei|+C^ 



7 / i - 1 



7-1' 



p+2 



T{a + l) f] 
T [a — p) \m 



7(a-p-l) 



^ ( L 

m\m 



7-1' 



P+2 . . ^ 7(a-p-l) 



As |ei| is less than some constant multiplied by the 'wth term' in the sum: 
C 



E„, 



7(a + l)-(p+2) 



(5) 
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Simplification of (||) (see below) leads to the result: 

Theorem 2.2 (Convergence of the h Method with Algebraic Grading) 

Consider the approximation of \x\°'dx, —1 < a < 1, using an h method based 
on a quadrature rule of degree p (an odd positive integer), on an algebraic mesh on 
a total of m ^ 2 intervals, with mesh parameter 7^1. For some constant C, the 
error Em satisfies: 

E,n S$ Cto"^. 

Here z is: 

^^f7(a + l) 1<7< (p+l)/(a + l) 
[ p + 1 else. 

When 7 = (p + 1) / (a + 1), £;,„ C In (to) /mP+i . 



Theorem 2.1 Take and write Em as: 



The integral converges if: 7(0; + 1) — (p + 2) < —1, and in this case, it converges 
to the constant [{p + 1) —7(0; + 1)] ^, independent of m. Absorbing this into the 
main constant, the result for 1 ^ j < (p + 1) / (a + 1) is created. 
In the case 7 = (p + 1) / (a + 1), is 



C C r 1 , C , , , 

Em < — - < —rr^ / -dx ^ —-3- In to) . 



Thus again, the error is bounded by a constant, this time, dependent on to. Ob- 
serve thatp+ 1 is immediately able to be replaced with 7 (a + 1), demonstrating the 
continuity of the formulae. 
Lastly, bound Em, as: 

C ™ / j \ •7(a+l)-(P+2) (J rl 



j=i 



l/m 



As 7 (a + 1) — (p + 2) > —1, the integral converges to [7 (a + 1) — (p + 1)] ^ , again, 
another constant independent of m, which is absorbed into C. The second result is 
thus achieved, by observing that to^'p^^) < m~^P~^^\ 

The moral of this is that for a particular choice of a and p, there is an ideal choice 
of 7, that is 7* = (p + 1) / (a + 1), beyond which the order of the error will not 
decrease. (Choosing 7 greater than this may reduce C.) Varying 7 makes no differ- 
ence to computational expense.^ Observe that 7* becomes unbounded as a — > —1^. 

^Choosing [7*] may be cheaper than nonintegral choices of 7*. 
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This is not surprising, as the integral itself becomes unbounded. A similar result to 



Theorem 2.2 exists for a geometric mesh. 

In summary, using an algebraic mesh of m subintervals, and an h method with 
a degree p quadrature rule on each mesh interval, can achieve, for integrands, 
an O (m~*-^^^^) convergence rate. This is a significant improvement on the equiva- 
lent result for a uniform mesh, which is O (m~'-"+^''). The exponential convergence 
rate for smooth integrands is not achieved, but can be, using an h-p method. The 
methods may be extended to integrands of the form \x\" f (x), for smooth functions 
/• 

2.6.2 Error Analysis for the h-p Method 

This section discusses the expected order of the error for the h-p method, for al- 
gebraic or geometric meshes using integrands with end point singularities. Again, 
consider the integral For the h method, Pj, the degree of the quadrature 

rule on the jth interval, was constant. For the h-p method, it is a function of j. A 
low degree rule is used on the interval adjacent to the end point singularity; and 
higher degree rules are used on intervals away from it. A simple choice is to use 
a rule on 2 points on the first interval, and increase the number of points linearly 
with j. That is, using Gaufi-Lobatto rules, where pj = 2nj — 3; choosing nj = .? + 1 
gives Pj -I. 

Recall the error result from (||), where Cj — C (pj) is given by (||): 

/ (x) dx - Ip^ = = CjhP+^f^P+^'> iO, e e , X,] . 

Use Stirling's formula for large x to approximate the factorials in Cj, as {x — 1)! = 
Tix): 

2tt rx\^ { ^ 1 



^^^^-\l-[-e) 12. 



Applying the first term of this to gives an asymptotic bound for large j: 

^ (j + i)(j + 2) [{j-m" 



Thus: 



^ = eV2^ (6) 

Em cannot be bounded directly using this approximation for Cj, as it is not a 
proper bound. Heuristically,^ for either a geometric or an algebraic mesh, the rapid 
convergence of Cj to means that Em is expected to be dominated by ei . Consider 
a mesh on m intervals, and an associated composite quadrature rule on a total of 
= TO (m + 1) /2 + 1 points. For a geometric mesh hi = cr™~^, and for an algebraic 
mesh hi = ra^^ . Using (^), this gives: 



This is brought out by experiment. 
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The error for an algebraic mesh is polynomial, whilst the error for a geometric 
mesh is exponential. If the errors with increasing N « rr? jl are plotted for a 
geometric mesh, an error of the form Em ^ Cp^ is observed (see Figure ||), for 
some p < 1. 

These rough results prompt more rigorous examination of Recall: 

For a geometric mesh, Xj = cr™^^, so hj < (1 — cr) a'^~K Using / (x) ~ gives: 
f{n),^^_ r(a + l) 



r (a + 1 - n) ' 



The Stirling asymptotic approximation for Cj can be converted to a genuine bound 
by observing for x € N: 

Recall from (||), that with fj, ~ for some C independent of 



IQI ^ ^' 



x: 

3 



pj+5/2 ■ 

A bound on the error for the jth interval is now:^ 



Simplification of this leads to a bound on Em- Observe that (1 — cr)^"'''''^ < 1, so: 

I I ^(r»-.7) (2.7+3) r(a+l) (,„_,-_i)(n_27-2) 

l^^l < ^j2j+5/2'^ r(a-2j-l) 
Absorb F (a + 1) into C, and use the Stirling approximation for F (q — 2j — 1): 
1 ^ /a-2j-l/ " ^"-^^-^ 



F(a-2j-l) V 27r V"-2j~l, 

Rearranging the exponent of cr, and absorbing the term e"^^ into C: 

< (7 '"^ ^(m-.7)(a+l)-a+2.7+2 . / « - 2j - 1 ^ " \ " 2j 1 



^ ^Jcr'"("+i)(jj(i-")e^2j 



j2j+5/2(^_2j_l)"-2j-l/2- 

Expand — e^/2^, and observe that {a — 2j — 1)" ^^'^ > 1: 

^ '^j2j+5/226j(„_2j-l)"'^' ^ 26jj5/2 

As I (a - 2j - 1) /j\ < 2, then: 

l^^l ^ 26jj5/2 24ij5/2 ■ 



^Strictly speaking, this oiily applies for j = 2 : m, as the result for Cj only holds for j = 2 : m. 
Application of the result in (m allows |ei| to be bounded by a constant multiple of this Ci. 
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Combine these, to bound: 



^ El 



24jj5/2 



^ 9 



i(l-a) 



^ 7-5/2 



Observing that 2^-' j^/^ > 1, this can be written as:^ 

m 

Using the result for the sum of a geometric progression: 



- 1 



1 



^(l-a)(m + l) _ 
fjl-" - 1 



Thus: S„ < Cfj^^^+iV^i-^^^^+i) < Ccr^". As w mV2, using p = cr^v^^ ^j^jg 
simphfies to: 

E„, < Cp^. 

This constant C is actuahy a function of a and cr, and the fact that the degrees of 



the quadrature rules are linearly graded. This proves Theorem 2.3 



Theorem 2.3 (Convergence of the h-p Method with Geometric Grading) 

Consider the approximation of Ixl^dx, —1 < a < 1, using an h-p method based 
on a geometric mesh on m intervals with parameter a, and using a Gaufl-Lobatto 
quadrature rule on j + 1 points on interval j . For some constants C , and p < 1, the 
error E„i satisfies: 

E„, < Ca^"". 



Babuska et al. [T2| |l^ describe an approximation theory for the h, p and h-p 
methods for the finite element method. This material may be able to be simplified 
and adapted (as the theory for integration should be easier than for approximation) , 
and also used in error analysis. 

2.6.3 Experimental Results for a Real Integral 

The example y/xdx is used to demonstrate the above convergence results. MAT- 
LAB code used (hpmeth.ni and funchp.m), is contained in Appendix and error 
results are presented in Figures |l|and| 

Figure |l| compares convergence rates for various choices of the algebraic grading 
parameter 7, whilst holding constant the number of points in the quadrature rule 
used on each interval; that is p = 6 is fixed. Figure ^ shows similar data, but varies p. 
(Here, 7 is allowed to vary, and is chosen to be equal to p for convenience.) Both plots 
are shown compared with the corresponding h-p result, using a geometric grading, 
and a = 0.15. The linear increase in number of points used in the quadrature rule 
means that Uj = j + 1, j = 1 : m, so A^ = m (m + 1) /2 + 1 (as rules are closed, but 

ends are not), and thus m cx Va". The h-p result demonstrates O (p^^) behaviour. 
The h-p data does not show as a straight line, but this can be seen in Figure ^, 
where it is plotted versus ^pN . 

*This throws away a lot of information! 



12 








loglO(points) 

Figure 1: Errors for h and h-p methods applied to ^/xdx, for various choices 
of 7 = 5. The slopes of the lines are approximately —87/2. As the error tends to 
machine precision (e w 10~^^), the convergence results lose their regularity. 

















P = a 


















Xn = 9 








h-p method 





0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 

loglO(points) 

Figure 2: Errors for h and h-p methods applied to ^/xdx, for constant p. The 
slopes of the lines are approximately — (p + 1). 
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2.6.4 Extension to Complex Contour Integrals 

The method is easily extended to complex contour integrals. Good test problems 
have closed contours, and integrals which can be directly evaluated using Cauchy's 
integral formula. To demonstrate this, consider / (z) = e~'^'^ I ^ \J z — 1, integrated 
around the unit circle. The integrand has a simple pole at z = 0, and a derivative 
singularity at z = 1. The resulting integral is: 

!{z)dz = e"/^. 

MATLAB code uscd (cint .m and f uncci .m) is contained in Appendix |a|. Error results 
for an h-p method using GauB-Lobatto quadrature rules are presented in Figure 
^. The mesh is geometrically graded, with parameter a — 0.15. For a segment of 
a closed contour, with a corner at either end, let D be chosen as the number of 
mesh intervals between each corner and a wide, central interval, so m = 2D + 1 is 
the number of mesh intervals over that segment. Here, as the contour is the unit 
circle, 2 artificial corners are placed, and D is varied from 8 to 15. The grading of 
the quadrature rules is similar to that used in 5:2.6.3 - the number of points used 



in the quadrature rule increases linearly with the number of mesh intervals from 
the nearest corner, starting at 2 on the interval nearest the corner, and finishing at 
Z? + 2 on the central interval. 



-10 



-13 - 



-14- 




12 14 16 18 20 22 24 

sqrt(N) 

Figure 3: Error results for the complex contour integral. 



Convergence is plotted for the logarithm of the error with VN. Observe that 
the plot is linear, that is, the error is O (p^^) • These superb results show that 
the method is excellent for the numerical approximation of closed complex con- 
tour integrals. This success motivates the use of the h-p method in the CBIEM, 
where quadrature rules for complex contour integrals are required in the numerical 
approximation of the solution to an integral equation. 
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2.7 Summary - Advantages of h-p Methods 

This section has discussed three important aspects of the numerical approximation 
of integrals with end point singularities: 

1. An h-p quadrature method is superior to other methods. 

2. A geometrically graded mesh is superior to other choices of grading (maybe 
only marginally better than an algebraic one). 

3. The appropriate family of quadrature rules to use is the Gaui3-Lobatto, as 
they are closed, and of maximal degree for the number of quadrature points 
used. 

The quadrature rule used in g3 is chosen in this manner. 
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3 The Complex Boundary Integral Equation Method 



3.1 Origins and Description 

The CBIEM is a technique which numericaUy approximates the solution of the 
Dirichlet problem.01t reformulates the solution of the Dirichlet problem as the real 
part of a function which can be found as the solution of a complex boundary integral 
equation. The solution of a discretised version of this integral equation is then found 
using a collocation technique. Finally, a discretisation of Cauchy's integral formula 
is used to approximate the solution to the original problem at interior points, based 
on the approximate boundary data. 

It is related to the 'Complex Variable Boundary Element Method' jisj, which 
is a Galerkin version of the same technique, using 'hat' functions as a basis. (The 
collocation method creates an approximation to the boundary data by interpolating 
from known data, whilst the Galerkin constructs an approximation in terms of a 
series of basis functions defined on segments of the boundary.) As originally stated, 
the CVBEM only works on polygonal domains,^^ whilst the CBIEM is more general 
in that it also works on non-polygonal domains. 

The Dirichlet problem considered is on an open, finite, simply connected and 
non-empty region f2 C R^. fl is bounded by F, a piecewise continuous, anticlock- 
wise oriented contour. F has a finite number of corners, at which its derivative is 
discontinuous. The Dirichlet problem is: 

Given boundary data /, find [/ : M subject to the conditions: 

V2;7(x) = 0, xef^, c/(x) = /(x), xeF. 

Thus, the problem is to find the solution to Laplace's equation over a region, 
given functional data around its perimeter. This has many applications in the so- 
lution of potential problems, such as electrostatics and fluid flow. The value of U 
at points interior to Q is determined by the boundary data being 'diffused' from 
the boundary inwards, according to the Laplacian operator. It turns out that the 
problem has a unique solution for all cases of /. In all but the most trivial cases, 
this solution is not expressible in closed form, and a numerical approximation is re- 
quired. With sufficient (possibly enormous) computational effort, an approximation 
to any degree of accuracy can usually be obtained. 

Problems with 'corner singularities' are of particular interest. In these problems, 
U is differentiable in the interior, but VC/ becomes unbounded as the corner is 
approached. It is known that this behaviour is typical of solutions to the Dirichlet 
problem on domains with corners. Even if the boundary data is smooth, VC/ still 
becomes singular near the corner. Numerical methods must be able to produce good 
approximations to U, in spite of the corner singularities. 

'Interior' methods, such as finite difference and finite element methods, become 
computationally expensive when applied to problems with corner singularities, and 
boundary integral methods are more appropriate. Interior methods require a finely 
discretised two dimensional mesh in the region of the corner, which greatly increases 
the size of the associated linear system. In contrast, a boundary integral method 
has only to discretise its mesh in one dimension, that of arc length on the boundary, 
and is expected to be much cheaper. 



^The space containing the functions approximating the solution of the Dirichlet problem is a 
Sobolev space, which is a generalisation of the Banach space of continuous functions to include 
functions with 'weak derivatives'. 
^"^The CVBEM has also been generalised to doubly connected domains [|lq|. 
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The usual boundary integral methods based on Green's functions lead to a kernel 
with a logarithmic singularity, even on a smooth domain. This is tedious to program, 
and computationally inefficient if high order methods are used. If the CBIEM is used 
with singularity subtraction, the integrands are smooth and can be done simply and 
accurately by direct quadrature. 

The problems caused by the corners and corner singularities are dealt with using 
h-p quadrature methods, and would be difficult to implement with other types of 
integral equations . 
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3.2 Development of the CBIEM 

The solution to the Dirichlet problem, U, is harmonic, as it satisfies Laplace's equa- 
tion in il. Identify x S with z G C Now U can be thought of as the real 
component of an analytic function W (z) ~ U (z) + iV (z), where V is uniquely 
determined to within a constant. V can be made unique by requiring V (^o) = for 



some Co € r (see § 3.2.5 ). For all z G T, U (z) = f (z) is immediately known. The 
CBIEM first approximates V (z) on F, and then uses Cauchy's integral formula to 
approximate W, and hence U, at points within CI. 

3.2.1 Cauchy's Integral Formula 

For an analytic function on a bounded domain fl, Cauchy's integral formula is: 
rxr (A\ fOz^fiur exterior 

—^dC = TTiW iz)x I 1 zeT boundary (8) 

^ ^ ^ \^ 2 z G interior. 

When z is on the boundary]"] the integral is a Hilbert transform. The kernel is 
singular, and the result must be interpreted as a Cauchy principal value integral 
page 39]. The CBIEM requires approximation of the Cauchy integrals by quadra- 
ture. In j |3.2.3 , this is used to set up a linear system for the approximation of V on 



r. After this has been done, in §3.3 it is used to compute an approximation to W 
(and hence U) in the interior of il. 

3.2.2 The Complex Boundary Integral Equation 

To derive the integral equation underlying the CBIEM, observe that letting W (() 
1 for the case z G F in (^) yields: 

dC = TTi, z E F. 

r ^ 

Multiplying this by the constant W (z) gives: 

W (z) i -r^dC ^ <j> ^^dC = T^iW (z) . 
JtQ- z Jr C- z ' 

Equating the niW (z) with that in (^) gives: 
W{C)-W{z) 



c- 



dC = 0, z e F. (9) 



(^ is called the 'Complex Boundary Integral Equation'. This derivation is par- 
allel to that involved in singularity subtraction page 184] and ]Q]. The integrand 
is analytic, and as C ^ z, it converges to W (z). An analytic function W — U + iV, 
which has as its real component the solution to the Dirichlet problem, will satisfy 
(^ . The converse is also true - a function W that satisfies (^ will have a real compo- 
nent U that satisfies the Dirichlet problem. It is hoped that a function that satisfies 
a discretisation of (^ will have as its real part the solution to a discretisation of 
the Dirichlet problem. 



^^This result is a simplification. If ^ is at a corner, replace 1 with a/ir, where a is the interior 
angle subtended by the corner (else a = tt). This result requires that collocation points are not 
placed at corners, to avoid unwanted complexities in the imn lemep tation. Fortunately, this is 



already overcome by the use of node points at the corners (see §3.2.3) 
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3.2.3 Discretisation of the CBIE 

The CBIEM requires the numerical approximation ol the (Cauchy) integral in (||), 
and this is achieved using quadrature. In particular, given the possibly singular 
nature of W at corners, §^ motivates the use of geometrically graded h-p quadrature, 
because it is of high order for such integrands. The approximation will be referred 
to as a discretisation. Nomenclature used in the following discussion is shown in 
Figure ^. (The distinction between mesh and node (quadrature) points is made in 
§ |S| .) 

segment 1 




non-smooth contour F = 6Q 



Figure 4: CBIEM nomenclature. Corner, mesh, node and collocation points on 
contour F, about a region O. The lengths marked on segment 1 are the positions 
of geometric mesh points, in terms of a unit arc length on that segment. The node 
points correspond to a closed Newton-Cotes rule on each mesh interval. 

Consider the contour integral of an arbitrary integrand g (C) around F: 

Parameterise F using 7 : [0,1] 1-^ C, such that Co = 7 (0) = 7(1) = (n, with 
argument t increasing in an anticlockwise direction around F. The contour integral 
is i page 168] :|3 



^^This requires that 7 is continuous, and that 7]^^ ^ j j is continuously differentiable for a finite 
partition = to < ii < . . . < tn = 1- 
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Approximate the integral using an h-p quadrature rule {tj , Wj } with A'^ node 
points, defining (^j — 7 {tj) and jj — — {tj): 

£ 9 (0 dC = l\ (7 (0) |J (t) dt^J2g {Q)i,w,. 

This formula can be applied to the Cauchy integrals. For some fixed z € T, consider 
the integrand 5(C) = {W {Q - W {z)) / {( - z). Let Wj = W {Q), and redefine 
Wj <J= 'jjWj to absorb 7^. The Cauchy integral is: 



W{C) -W{z) 



N 



Wj -W{z)^ 
Cj-z 



z e r. 



(10) 



Approximation of the solution to the CBIE requires approximation of the Cauchy 
integrals without using W {z). Instead, W {z) k, W {z) is constructed from values 
of W at the quadrature points. (|o|) is discretised into a linear system of order A^. 

Begin by choosing a set of A^ different values of z from around the boundary. 
These points are called the collocation points. It is known from analysis in the 
case of uniform meshes that collocation points must not lie on node points [ pO[ . 
The natural choice is to take as the collocation points the A^ midpoints (in the 
sense of arc length) between the A^ node points.]^ Let tfc_i/2 = {tk-i +tk) /2, 
Cfc-1/2 = 7 {tk-1/2) and Wk-1/2 = W (C/c-1/2), for fc = 1 : A^. Interpolation from 
known values of U at the node and collocation points is used with the CBIE to 
approximate the Wk^i/2- 

Define two complex A^-vectors of W at the node and collocation points: 



W = 



Wi 



Wi 



N 



w - 



w, 



N-1/2 



Also define the real A^-vectors U = (W), U' = (W) and V = Q (W). From 
( p^ , an order A^ linear system for the components of W and W is determined: 

N 



E 



Wj - Wk^y2 

Q - C 



*:--l/2 



Wj = 0, 



fc = 1 : A^. 



(11) 



In summary, discretisation of Cauchy's integral formula leads to a linear system, 
the solution to which is an approximation to W at the A^ node points on F. This 
approximation can be used to approximate W, and hence U, at points within fl. 



'^^The choice Ck—1/2 = (Cfe— 1 + Cfe) /2 is explicitly not used, as this assumes the contour is linear 
between points parameterised by and t^. 
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3.2.4 Linear Interpolation of W at the Collocation Points 

If the Wj were known, ( |Tl| ) could be directly used to interpolate the VFfe_i/2- How- 
ever, although Uj and ?7fc_i/2 are known explicitly, Vj and Vfc_i/2 are not. The 
CBIEM implicitly approximates the Vfc_i/2 by interpolation from the as yet un- 
determined Vj at points near Cfc-i/2- That is, if a rule on O points is being used, 
choose O terms from the sequence: 

Cfc-3, Ck-2, Cfc-i; Ck, Ck+1, Ck+2, ■■■■ 



The discretisation of the CBIE in (|ll[), coupled with the 2N knowns Uj and Uk-1/2, 
and the interpolation for the Vk-1/2, allows the approximation of the N unknowns 
Vj. (The Vk-i/2 are not explicitly required to be calculated.) 

The simplest interpolation for the 14_i/2 is a linear one, between Cfe-i and Ck, 
that is use: 



14-1/2 « (14-1 + 14) /2. 



(12) 



More sophisticated interpolations to the Vfe_i/2 could involve higher degree poly- 
nomials, splines, or trigonometric polynomials. Initially, the linear choice will be 
used to illustrate the process. In § ^.5.1 , the method is extended to higher degree 
polynomials. For a particular problem (fi, /), there is an optimal choice of degree 
for the interpolation, as errors incurred by interpolation increase with the degree, 
and eventually become of greater magnitude than those due to discretisation. 

The Wj are found by solving for their imaginary parts Vj , and combining these 
with the knowns Uj. The linear system is set up as follows. Substituting Wj — 
Uj + iVj and M4_ 



1/2 

knowns and unknowns: 



Uk-1/2 + «14-i/2 into (|ll|), using (jlj), and collecting 



^ '^Vk-i + IVk 



V 



Cj ^ Cfe-1/2 



N jj 
.y^ L/fe-1/2 

j^l Cj ^ Cfc-1/2 



k = l: N. 



(13) 



This is a system of N equations for the N unknowns Vj, with RHS determined by 
the knowns Uj and Uk-1/2 (and of course the associated Q and Cfc-i/2)- 



3.2.5 Solution of the Collocation Equations 

In order to write (|l3|) as a linear system, consider the LHS of its fcth equation: 



N 



N 



j^i - Cfe-1/2 

Define a matrix A e C^""^: 



N 

E 



j~[ - Cfe-1/2 ^ Cj - Cfc-1/2 



W2 



WN 



Cl - Cl/2 

Cl - C3/2 



Wl 



C2 - Cl/2 
W2 

C2 - C3/2 



■W2 



Cn - Cl/2 

Cn - C3/2 



WN 



Cl ^ CAf-1/2 C2 — CAf-1/2 Cw ~ CjV-1/2 

Also define a set of N scalars Hk, for fc = 1 : TV (the row sums of ^4): 

N 

E7/I.- 



Cj - Cfc- 



1/2 



(14) 
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The first two terms of (|lj) are then ^Vk-iHk 
and hence Vq^Vn- Define B e C^^^: 



\VkHk. As r is closed, (o = C 



N 



Hi 
H2 



Hi 



H2 



H 



Hn-1 
Hn 



H 



N 



(15) 



Let C = B - ^, then the LHS of (y|) is CV. Let In G represent the real 
column vector with all components unity, and the operation diag (x) on A^- vector x 
create the diagonal matrix of order N with the diagonal entries being the respective 
components of x. Defining d = diag (Aljv) U' — AXJ, the system for V is: 

CV = id. (16) 



Attempting to directly solve the complex linear system in (16) fails, as V is overde- 
termined in two separate ways. Firstly, V is purely real, that is 3 (V) = 0. Par- 
titioning (^6|) into real and imaginary components yields two purely real linear 
systems, either one of which can be solved for what should be the same solution V. 
The system to be solved is: (C) V = -3 (d) or 3 (C) V = (d). By redefining 
C <;= 3?(C) and d —3 (d), the first choice for the solution of V is made. The 
linear system is thus: 

CV = d. (17) 

The second way that (|l^ ) is overdetermined is that V is known only to within 
a constant.0 So, if direct solution of (|l7|) i s attempted, singularity problems will 
occur, and the resultant V will be infinitep] To accommodate this, arbitrarily pi se t 
Vn = 0, and compute the rest of the components of V by subtracting rows in (|l7|). 
The result is a fully determined order A^ — 1 linear system. The A^ rows of (|l7| ) are: 

(CV)^. -d„ j-l:A^. 

Subtracting rows gives a system of A^ — 1 equations: 



Defining V* e K 
C*V* = d' 

Here: 

C* = 



(CV),-(CV),_i-d,-d,_i 



N- 



'J \ 'J 
as the first A^ 



J 



2 : N. 



1 entries of V (where the last entry is zero): 

(18) 



C2,l — C14 

Ca^i — C24 



C2, JV- 
C's, jv- 



1 — Ci^N- 
1 — C2.N- 



d* = 



Cjv,i 

d2- 

d3- 

(In — d 



- Cn^ 

di 
d2 



N~l 



Cn- 



1,JV-1 



b(JV-1)x(JV-1) 



Solution of the order N — 1 linear system in (|18|) yields the approximation to V, 
and hence W {W at the A^ node points). 



'^^A scalar multiple of Ijv- 

'^^Well, a numerical approximation to 00! 

^^For test problems, Vjv = V (Cn) is actually used. 
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3.3 Approximation of U at Interior Points 



The approximation W = U + iV on F is used to approximate W (and hence U) at 
interior points of Q, using Cauchy's integral formula for points within Q: 



W{z) 



1 



Discretising this gives the approximation, for z e fi: 

N 



W(z)^J-^ 



2ni ^ Q 



(19) 



Numerical problems occur using this simple approximation, as points z G near the 
boundary (where z — Q is small, for some <^j), generate very large terms in the sum. 
In fact, the integrand is nearly singular, so only poor acc uracy is expected. Instead, 
the technique of singularity subtraction (referenced in § 3.2.2| ) uses the result from 



1 ^ 

— T 



Wi Q)-W{z) 



Wj = 0. 



The integrand is now smooth, and good results can be expected from quadrature. 
This yields W {z) as a 'corrected' (noh: 



W{z) = 



N 

E 



0-^ 



N 



1 







Implementation of the CBIEM using this result is successful. The code supplied 
(see does not go beyond the stage of the computation of V on F, as it is known 
that the approximation to U in the interior is actually more accurate than the 
approximations to V on the boundary Experiment demonstrates this, and thus 
computation of U at interior points need not be further described. 
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3.4 Performance of the CBIEM on Model Problems 



Application of the CBIEM yields different quality results depending on the conti- 
nuity of the model problem and the contour. In the simplest case, both are smooth, 
there is no singularity, and standard quadrature gives good results. In fact, because 
of periodicity, even the trapezoidal rule on a uniform mesh gives very good results. 
There is no need to use h-p quadrature, but it will still work well. 

Now consider the case where F is smooth, but W has singularities. For exam- 
ple, if VF(7(s)) ~ the integrand of the CBIE §) can have behaviour s ' 
near the corner. However, using singularity subtraction and h-p quadrature, exper- 
iment demonstrates that both the discretisation error and the final error in V on 
the boundary are superb. (See also the example of complex contour integration in 
§||6j.) 

If W is smooth, but F is not (has corners), in general, the errors will be expected 
to increase with the sharpness of the corner. The most difficult cases are cusps or 
reentrant corners (e.g. the corner in a cardioid). Even for a model problem with a 
smooth solution U, a corner singularity in V (and hence W) will occur. 

Let r represent radial distance from a corner. It is known pages 257-259], 
that at a corner with interior angle (1 — x)'^t ^ singularity of the form r^/(^^x^ 
will be found. At worst, for a reentrant corner, x = ~lj so the form is r^/^. A 
good model problem is thus a contour with a corner where the true solution has 
local behaviour U ~ r^^^. This is obtained, for example, using W (z) = z^l"^. If a 
uniform mesh were used, the greatest component of the error in V will come from 
the intervals adjacent to the corner. The use of a geometrically graded mesh reduces 
this component to a level comparable with that of other mesh intervals. Errors will 
not be of the very high order that is expected for smooth contours, but should still 
be acceptable. 
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3.5 Higher Degree Interpolatory Polynomials 
3.5.1 Introduction 



To extend the technique described in §3.2.4, the hnear interpolation in <\12) is re- 



placed by a higher degree interpolation. The net result of this is to change the 
definition of the matrix B in (p^). Other possible techniques of improving the ac- 
curacy of the interpolation, such as splines, are not considered here, as they are 
difficult to implement. The principle involved is that increasing the order of the 
interpolatory polynomial should reduce the discretisation error, which is expected 
to be greatest on the largest intervals. 

The 'nearest' O node points on either side of Cfc-i/2 are used to yield an interpo- 
latory polynomial of degree O—l. In general, for points far from the nearest corner, 
this means to take the first O terms of the sequence Cfc, Cfc-ii Ck+i, Ck-2, ■ ■ ■■ Oth- 
erwise, the term 'nearest' is used loosely, as interpolation cannot continue around 
a corner. Where there are less than 0/2 node points between the collocation point 
and the nearest corner (including the node on the corner); instead O points from 
and including the corner are used. This results in an interpolatory polynomial that 
is expected to be least accurate at the collocation point adjacent to the corner. 
(A possible improvement in this schema is to organise the interpolation rules such 
that their order increases say, linearly, with node index away from the corner.) The 



constraint on O due to the mesh parameter D is described in §3.5.4. 

§§3.5.2 and 3.5.3| deal with technical implementation issues, and can be skipped 



without loss of continuity. 

3.5.2 Lagrange Form of the Interpolatory Polynomial 

The Lagrange form of the interpolatory polynomial is used. Given a set of values 
for the Vj, at positions Q — 'y{tj), the approximation to V at point Cfe-i/2 is: 

Here:[] 

A, w= n 



tj 



is a set of indices of the nearest node points, specific to the collocation point 
C/c-i/2- Specifically, where a degree O — l interpolatory polynomial is used on O 
points, and C is the index of the nearest corner to Ck-1/2' 

( {k - 0/2,..., k-l,k,...,k + 0/2-1} in general 
Fk= I {C,...,C + 0~l} iik-0/2<C 
\ {C -0 + l,...,C} iik + 0/2>C. 

Let Fk be the kth row of a table F, and let the above \j be the jth element of the 
kth row of another table, L, of the associated weights. The notations F{k,j) and 
L {k,j) are used to describe the set of O nodal indices and weights associated with 
the interpolation at point Cfc-i/2i where fc = 1 : A'^ and j = 1 : O. The structure 
of F becomes more complicated with increasing O and with increasing number of 
corners. Details of its construction are not provided here, but can be read from the 
program cbiem.m in Appendix [A.l| . 

'^^Warning: Replacing t,y with (^^ and tj with {^j in this formula cannot be done, as the contour 
segments are not necessarily straight. 
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For some integer D, on each segment of F (with a corner at each end), a mesh 
is constructed that has 13 — 1 internal points between each corner and a (wide) 
interval which spans the centre of the segment. This results in a total of 2D mesh 
points, and 2D — 1 mesh intervals (see Figure ^ . 

For example, if D = 3, then there are 4 interior and 2 corner mesh points on 
each segment, together with 4 extra interior points, for a total of 5* = (D + 1)^ + 
1 10. Further, if there are NC = 3 corners, then the linear system has order 
N = NC{S — 1) = 27. If O = 6 (quintic interpolation about the 'nearest' 6 points), 
then F E ^Nxo_ \Yhere divisions in the structure of F due to the corners are 
reflected by partitions, F isiPl 



F = 




Ft{l) 
Ft{2) + ones{Ft) 

Ft{NC) + {NC - l)^ * ones(Fi) 



For j = 1 : 5 - 1 and fc = 1 : O, in general Ft e N^'^"^)^'^ is: 
Ft{j,k)^ 

(Exception: F(l : 0/2, 1) x lo/2-) 



1 J < 0/2 

l+j-0/2 0/2<j<S-0/2 
l + S-O S'-0/2«;j 



3.5.3 Construction of B 

The matrix B is required in the construction of the linear system in (p^, and is 
the only thing that changes when O is varied. For the case of linear interpolatory 
polynomials {O — 2), a formula involving the terms is used to approximate the 
value of V at the collocation points. For larger O, this is replaced with a considerably 
more sophisticated formula. As O increases, the bandwidth oiB increases. Naturally, 
this new B simplifies to the earlier definition if O = 2 is used, but is obtained at 
greater computational expense. 

The crucial change is in the approximation to the value of Vk-l/2^ which in ( p^ ) 
is the term \Vk-i -\- \ Vk. This is now replaced with: 

o 

^^Liberal use of MATLAB notation is made, and tliere is a confusion between computer array 
element and mathematical matrix entry notation: F {j, k) = Fj f^. 
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The first O terms of the LHS of the fcth hne of ( p_3| ) (there is only one other term) 
are now: 



o 

E 



L {k,j) VF{k,j)Hk. 



Construction of the real order N matrix with (fc,j)th entry L {k, j)Vp(^k j) is re- 
quired. Multiplying each row by the corresponding converts this to B. 

Details of the construction of B are not provided here, but the illustrative exam- 
ple used in §3.5.2 is continued. Recall that D = 3, NC = 3 and O 
and N 



6, so 5' = 10, 



27. First construct Bt (a 'skewed' version of L), and then calculate B by 
multiplying Bt through by the Hk- Bt is constructed using a 'shift vector' G, where 
Gk is equal to the number of zeros to be put in front of row A: of L to make it row 
k of Bt. This G has a structure formed from a temporary Gt: 



Gt = 

G = 
Bt{k,:) = 



zeros (1,0/2) [1 : 5 - O - 1] (S-O)! 

T 



T 

0/2 



[Gt Gt + S-l ... Gt + NC{S -1)^ 

[zeros (l,G(fc)) i (fc, :) zeros (1, TV - G (fc) - 0/2)]^ k = 1 : N 
B <^ Bt{:,2:N) B {1 : 0/2, N) ^ Bt (1 : 0/2,l) . 



The overall structure of B is depicted in Figure ||, where x and • represent 
nonzero and zero entries respectively. 



X X X X 



X X X X 



X X X X 
X X X X 



X X X X 



X X X X 



Figure 5: Structure of B. 

The NG X NG submatrices within the structure are each of order S — 1. The 
kth row generally consists of O = 6 contiguous nonzero entries, starting at column 
Gk. These entries are the kth row of L, multiplied by Hk. That is, the string of 
O = 6 nonzero elements in row k represents: 



[HkL 



fc,i 



HkLk,o-i HkLk^o] 
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Additionally, in the first 0/2 rows, the first of these O entries is shifted to the 
A^th column, due to a 'wraparound' effect. In view of the banded structure of B, it 
appears that a method designed to exploit this structure would be appropriate in 
the solution of (|8|). Unfortunately A is neither banded, nor of a particularly simple 
structure, so this approach is nontrivial, and is a direction for further work. 

3.5.4 Limitations on the Degree of the Interpolatory Polynomial 

The discussion and results presented in §^ prompt the use of a geometrically graded 
mesh, with the number of points in the (closed) quadrature rule on each mesh 
interval linearly increasing with interval number from the corner, beginning with 2 
adjacent to the corner, and becoming D on the central (widest) interval. The use 
of a linear grading is found in the literature of the finite element method and the 
usual boundary element method [Q, ^] . Changing from linear to quadratic or higher 
degree may reduce the errors, but its implementation is beyond the scope of this 
report. 

Consider a segment of F divided into a mesh on 2 (_D + 1) points, including 
corners, with m = 2D + 1 intervals. Basic quadrature rules on rij — 2,5, D — 
1,D,D — 1,...,3, 2 points, are used over intervals j — 1 : m. After all of the 
common end points are considered, the total number of points in the final composite 
quadrature rule for that segment is S = {D + 1)^ + 1. Slightly different results would 
apply if the quadrature rules were open. Recall that this choice is rejected, as it 
leads to node points that avoid the singularity. 

S imposes a limit on O, the number of points used in the interpolation rule. As 
O is usually small (for reasons of computational efficiency, typically O ^ 6), this 
limitation is usually not significant. For example, ii D = 2, O ^ 8, and if D = 3, 
O ^ 16. 



28 



4 Implementation and Results 



4.1 Implementation 

The CBIEM is implemented as a set of functions in MATLABp] code, presented in 
Appendix^ (These functions appear in alphabetic order, interspersed with several 
functions referenced in §|^.) The main routine is cbiem.m. The parameterisation 
of the contour and its derivative are computed within cbiem.m, and it calls an 
auxiliary function (funccb.m) to compute the true solution for a test problem. The 
quadrature points for basic rules are obtained by calling the function gettw . m, which 
provides either (closed) Newton-Cotes points (using a routine internal to gettw. m), 
or calls another function, lobatto .m, which computes the points for a closed Gaufi- 
Lobatto rule. Two further functions arc used to create an h-p composite quadrature 
rule out of a set of basic quadrature rules (hprmesh.m), and compose a quadrature 
rule over a contour with several corners (rmesh.m). For testing cbiem.m over a large 
set of parameters (e.g. generating the data for Tables |^ to a driver function, 
testcb.m, is used. 

Within cbiem.m is a description of its input parameters. For test problems, 
where the true solution is known, it plots and calculates norms of V — V, and 
also calculates some discretisation errors. Explicit computation of the approximate 
solution within is not performed. Experiments with doing this demonstrate that 
the error results obtained are of the same order as those returned. 

^^MATLAB is an (interpreted) matrix computation package, and is a trademark of The Mathworks, 
Inc. 
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4.2 Experimental Results using a Teardrop Contour 
4.2.1 Description 

Although the CBIEM code is generahsed to the situation of multiple comers, good 
experimental contours have only one corner, to facilitate isolation of the sources 
of error. This section describes numerical results for the CBIEM, using a teardrop 
contour,^ depicted in Figure |^. The contour is parametrically given by: 

7 (t) = 2 sin (TTt) + i sin (27rt) t G [0, 1] . 

This is the same contour as that used in It has a right angle corner at the 
origin, which facilitates the use of test problems W {z) — z". For a € (0, 1), there 
is a discontinuity in the derivative of the true solution at the origin, which becomes 
more pathological as a ^ 0. 



0,8 
0.6 
0.4 
0.2 

-0.2 
-0.4 
-0.6 
-0.8 




Figure 6: Teardrop contour used in the CBIEM experiments. 

Error results are presented using an unweighted vector 2 norm: 

l|V-V||,= (K:-V-)']'^'. 

An appropriately weighted discretisation of the L2 norm might seem more appro- 
priate, but would effectively only present the norm over the central interval, as the 
widths of the end point intervals are very small. The use of an infinity norm is also 
appealing, but the 2 norm allows the user to experiment with interpolation formula 
gradings, to independently reduce the error over different regions of the contour (see 
4.2.6| ). Also, experimental data shows the behaviour of the infinity norm is very 



similar to that of the 2 norm. 

Tables |^ to |l0| present error results for the teardrop contour, for three different 
model problems: W {z) = z^, z^/^ and z^/^; various choices of two mesh grading 
parameters (cr and D); and choices of O, the number of points used by the interpo- 
latory polynomial in the collocation process. In each table N = {D + 1)'^ is the size 
of the linear system being solved, such that there are 2D + 1 mesh intervals between 



one corner and the next (see § |3.5.2|) . The nine tables cover three illustrative choices 
of the mesh parameter a for each of the three model problems. In each case, the 
results presented are for a choice of a close to the optimal a, and two nearby values 
of a that demonstrate the increase in the error in each direction. Results have been 
selected from a much larger data set. Within each table, the minimum error result 
is emboldened. 



-^Another important test contour is a cardioid with a reentrant corner. 
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4.2.2 Observations of 



This test function does not have a singularity, and the results are good. Despite 
the corner, the error reduces with increasing either D or O, until a point is reached 
where roundoff error, caused by excessive order in the interpolatory polynomial, 
begins to encroach. 





o 


N 


8 


10 


12 


14 


16 


18 


20 


22 


25 


4.6C-04 


5.6e-05 


6.GC-06 


l.Oc-06 


1.2C-06 


1.2C-06 


3.7C-06 


3.GC-06 


36 


1.4e-04 


1.3e-05 


l.le-06 


8.0e-08 


4.2e-08 


2.8e-08 


7.8e-07 


4.1e-07 


49 


7.1e-05 


5.7e-06 


3.9e-07 


2.5e-08 


5.6e-09 


3.5e-10 


5.8e-08 


2.9e-05 


64 


3.1C-05 


2.1C-06 


1.2e-07 


6.5C-09 


1.5C-09 


9.0C-09 


2.6C-07 


5.2C-02 


81 


1.6e-05 


l.Oe-06 


5.5e-08 


2.5e-09 


3.6e-10 


1.2e-08 


3.6e-05 


2.0e-02 


100 


8.6e-06 


4.9e-07 


2.3e-08 


9.6e-10 


8.4e-10 


l.Oe-08 


4.2e-05 


1.6e-01 



Table 2: ct = 0.20, [/ (z) = 3? {z^) . 





o 


N 


8 


10 


12 


14 


16 


18 


20 


22 


25 


5.9e-05 


6.6e-06 


9.0e-06 


9.3e-06 


9.3e-06 


9.4e-06 


9.3e-06 


9.7e-06 


36 


1.5e-05 


9.3C-07 


1.8e-07 


2.2C-07 


2.3e-07 


2.3C-07 


2.4e-07 


1.4e-06 


49 


6.2e-06 


3.1e-07 


1.2e-08 


5.0e-09 


5.5e-09 


5.5e-09 


6.9e-09 


6.7e-08 


64 


2.6e-06 


l.Oe-07 


4.0e-09 


9.9e-ll 


1.2e-10 


1.3e-10 


1.5e-10 


4.7e-09 


81 


1.3C-06 


4.5e-08 


1.4e-09 


4.1C-11 


2.3e-12 


3.1e-12 


4.6C-12 


8.1e-09 


100 


6.7e-07 


2.0e-08 


5.7e-10 


1.4e-ll 


3.8e-13 


5.2e-13 


2.5e-ll 


4.6e-08 



Table 3: a = 0.28, U{z) = ^ {z'^) . 








N 


8 


10 


12 


14 


16 


18 


20 


22 


25 


1.3e-04 


1.4e-04 


1.4e-04 


1.4e-04 


1.4e-04 


1.4e-04 


1.4e-04 


1.4e-04 


36 


6.3e-06 


6.7e-06 


6.8e-06 


6.8e-06 


6.8e-06 


6.8e-06 


6.8e-06 


6.8e-06 


49 


3.4C-07 


3.1e-07 


3.2C-07 


3.2e-07 


3.2e-07 


3.2C-07 


3.2e-07 


3.2e-07 


64 


l.le-07 


1.3e-08 


1.5e-08 


1.5e-08 


1.5e-08 


1.5e-08 


1.5e-08 


1.5e-08 


81 


5.4e-08 


7.9e-10 


6.8e-10 


6.9e-10 


6.9e-10 


6.9e-10 


6.9e-10 


6.9e-10 


100 


2.6e-08 


3.7e-10 


2.8e-ll 


3.1e-ll 


3.1e-ll 


3.1e-ll 


3.1e-ll 


3.1e-ll 



Table 4: ct = 0.35, U{z) = ^ {z'^) . 
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4.2.3 Observations of z^/^ 



This case has a corner singularity, and represents the 'worst' that singularities get 
in practice (that is, for z" singularities, in practice a ^ 1/2). Although the error 
decreases with increasing D or O, it does so more slowly than for z^, and is orders 
of magnitude larger. As for z^, there comes a point where increasing O causes the 
error to increase, and indeed grow exponentially. The results for V — V for the case 
cr = 0.10,r' = 9,O = 6 are plotted in Figure 0. The abscissae are plotted uniformly, 
for if they were plotted versus parameter t, the geometric grading would bunch up 
most of the results at the ends (corner of teardrop). Observe that these results are, 
as expected, antisymmetric. 





O 


N 


2 


4 


6 


8 


10 


16 


2.2e-02 


l.le-02 


2.6e-02 


7.3e-01 


2.3e+00 


25 


1.4e-02 


6.6e-03 


8.5e-03 


1.8e-01 


4.5e-l-00 


36 


l.Oe-02 


2.2e-03 


2.1e-03 


4.3e-02 


5.3e+00 


49 


8.1e-03 


1.3e-03 


9.3e-04 


9.7e-03 


l.Sc+OO 


64 


6.3e-03 


6.8e-04 


3.6e-04 


2.2e-03 


4.5e-01 


81 


4.9e-03 


4.5e-04 


1.9e-04 


5.5e-04 


l.Oe-01 


100 


3.9e-03 


2.9e-04 


l.le-04 


1.7e-04 


2.2e-02 




Table 5: cr = 0.05, U (z) = 


3fi(zi/2). 






O 


N 


2 


4 


6 


8 


10 


16 


2.0e-02 


1.4e-02 


1.3C-02 


4.4C-02 


2.7e-01 


25 


1.2e-02 


5.2e-03 


5.1e-03 


1.4e-02 


1.4e-01 


36 


8.3e-03 


1.6e-03 


1.5e-03 


4.5e-03 


4.9e-02 


49 


5.9e-03 


6.2e-04 


5.2e-04 


1.4e-03 


1.5e-02 


64 


4.5e-03 


2.5e-04 


1.6e-04 


4.5e-04 


4.9e-03 


81 


3.4e-03 


1.3e-04 


6.0e-05 


1.4e-04 


1.5e-03 


100 


2.7e-03 


8.4e-05 


2.4e-05 


4.6e-05 


4.9e-04 




Table 6: cr = 0.10, C/ (z) = 


3fi(zi/2). 






O 


N 


2 


4 


6 


8 


10 


16 


2.4e-02 


2.3e-02 


2.2e-02 


2.8e-02 


3.0e-02 


25 


1.2e-02 


9.7e-03 


9.3e-03 


l.le-02 


1.3e-02 


36 


7.0e-03 


3.7e-03 


3.5e-03 


4.3e-03 


5.4e-03 


49 


4.7e-03 


1.4e-03 


1.3e-03 


1.6e-03 


2.1e-03 


64 


3.4e-03 


5.7e-04 


5.3e-04 


6.5e-04 


8.2e-04 


81 


2.6e-03 


2.3e-04 


2.0e-04 


2.5e-04 


3.1e-04 


100 


2.0e-03 


9.8e-05 


8.0e-05 


9.8e-05 


1.2C-04 



Table 7: cr = 0.15, C/ (z) = 3? (z^/^) . 
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4.2.4 Observations of z^l^ 



A z^/'* singularity is beyond the range of singularities expected for smooth test 
functions. The errors are worse again than for z^/^, and they do not decrease as 
fast with increasing D or O. In fact, when the test problem is this pathological, 
the Dirichlet problem is fast becoming a boundary layer problem, which should be 
dealt with using more specialised methods. 





o 


N 


2 


4 


6 


8 


16 


4.6e-02 


4.4e-02 


8.1e-01 


l.Oe+01 


25 


3.2e-02 


3.3e-02 


2.8e-01 


3.2e+01 


36 


1.7e-02 


1.9e-02 


l.Oe-01 


4.9e+01 


49 


l.Oe-02 


1.2e-02 


4.1e-02 


3.3e+01 


64 


5.9e-03 


6.9e-03 


1.6e-02 


1.3e+01 


81 


3.9e-03 


4.2e-03 


7.5e-03 


1.8e+01 


100 


2.9e-03 


2.5e-03 


3.6e-02 


2.8e+02 


Table 8: a 


= 0.02, U{z) = ^{z 






O 


N 


2 


4 


6 


8 


16 


3.6e-02 


3.8e-02 


5.8e-02 


1.4e+00 


25 


2.0e-02 


2.1e-02 


3.0e-02 


7.5e-01 


36 


l.Oe-02 


9.8e-03 


1.4e-02 


3.6e-01 


49 


5.4e-03 


4.9e-03 


7.1e-03 


1.7e-01 


64 


3.2e-03 


2.3e-03 


3.3e-03 


8.3e-02 


81 


2.1e-03 


l.le-03 


1.6e-03 


4.0e-02 


100 


1.6e-03 


5.6e-04 


7.8e-04 


2.0e-02 


Table 9: a 


= 0.05, U{z) = ^{z'^/^). 




O 


N 


2 


4 


6 


8 


16 


5.8e-02 


6.0e-02 


5.8e-02 


6.4e-02 


25 


3.4e-02 


3.5e-02 


3.4e-02 


3.7e-02 


36 


1.9e-02 


1.9e-02 


1.9e-02 


2.0e-02 


49 


l.le-02 


l.le-02 


l.Oe-02 


l.le-02 


64 


6.2e-03 


6.2e-03 


6.0e-03 


6.5e-03 


81 


3.6e-03 


3.5e-03 


3.3e-03 


3.7e-03 


100 


2.1e-03 


1.9e-03 


1.8e-03 


2.0e-03 



Table 10: a = 0.10, U{z) = ^ (z^/^). 
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4.2.5 The Black Art of Choosing a 



The minimum error results obtained for each test problem are plotted versus a 
in Figure ^, and demonstrate that there is an optimal choice of cr, which varies 
significantly with the test problem. The use of the CBIEM in applications, where 
the true solution is not known in advance, could falter on the setting of a. If the 
computational cost is to be minimised, then it is important to find the optimal 
cr, however, it may be expensive to try many a until the optimal one is found. 
The literature does not justify a choice of cr, but merely states it, e.g. |Q uses 
a = 0.15 for a particular (finite element) application. The optimal choice of a for 
the paradigm test problem z^/^ is cr = 0.10. As z^/^ is the worst singularity expected 



in practice (see §3.4), this should be a good guide as a starting guess for any problem 
with an unknown solution. 



















y*-*-*^ z'^{l/2) 
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-4 - 



-6 - 



-10- 



-14 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 

Sigma 



Figure 8: Variation of minimum error with a . 



4.2.6 Improvements in the Technique 

Consider Table ^, where the best error result is obtained using O = 6. The error 
may be able to be reduced by grading the order of the interpolatory polynomial 
over the mesh intervals. Near the corner, the use of high order interpolation may 
actually increase the component of the error, although this may be appropriate 
far away from the corner. It may be ideal to grade the order of the interpolatory 
polynomial from O = 2 near the corner, to O = 6 (or greater) farthest from the 
corner. 
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Direct implementation of this result requires extensive modification to the matrix 
B used by the CBIEMQ (see § ^.5.l| ), and is beyond the scope of this report. Another 
way of achieving the same effect is to calculate B matrices B2, B4, . . . , Biq for 
= 2,4,.... 16, then construct a new B from appropriate rows of them, and insert 
this new B at the relevant point in the CBIEM. 



4 - 
2 - 

0- 

-2 - 
-4 - 
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30 



40 
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Figure 9: The best error results for z^/^, D = 9, cr = 0.10 and a graded interpolation 
rule. c.f. Figure]^. 



However, intuition is misleading here. The minimum error in Table ^ is 2.4158 x 
10~^, using a = 0.10 and D = 9. This corresponds to O = 6 on each mesh inter- 
val. Many experiments in variation of the order of the interpolation rule, hold- 
ing fixed a and D, find that the very best error result that can be obtained 
is 1.7444 X 10~^ (a 28% reduction), using a grading with interpolation rules of 
O — 12, 2, 2, 6, 6, 6, 10, 10, 10 over the 9 mesh intervals from the corner to the cen- 
tre. Surprisingly, the component of the error over the first interval decreases with 
increasing O. This is depicted in Figure ||, which shows that the error is uniformly 
distributed around the contour, except for the largest component, at the corner. 

It appears that what is happening is that the method has come up against 
a discretisation error barrier. For this problem, the discretisation error does not 
decrease particularly quickly, and is a maximum at the corner. 



^^These comments also refer to the associated F and L matrices. 
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5 Further Directions for Research 

This section enumerates various possibilities for future work on the CBIEM. 

1. There is the potential for error reduction using graded interpolation rules (see 
§4.2.6). Similarly, other choices for the grading of the quadrature rules may 
assist in error reduction, e.g. quadratic increase in degree of quadrature rule 
with node number from the corner, rather than linear as is presently used. 

2. A proper examination of the computational efficiency of the CBIEM is re- 
quired. This would involve setting up, say, a finite difference solution for the 
Dirichlet problem, and comparing flop counts required to obtain comparable 
accuracies. 

3. Analysis of the choice of optimal a is desirable. Currently, the method is 
hampered by this not being known in advance. 

4. Application of the technique to conformal mapping 0| may be worthwhile. 

5. It would be computationally efficient if solution of the linear system involving 
the matrix C = B — A could exploit the banded structure of B (see §3.5.3). 

6. An alternative collocation technique is possible |2^. Given nodes with 
weights Wj, and collocation points Cj+i/2 with weights 

. N N 

Approximate the unknowns Vj = ^(Ci) by collocating at Cfe-i/2j and the 
unknowns Vfc_i/2 = V (C/c-1/2) by collocating at C,j. This gives an order 2A^ 
system for the 2A^ unknowns Vj and Vfc_i/2, but avoids interpolation. 

■^1 - U-1/2 



N 

k=i '''=-1/2 - 4j 



The technique appears to be computationally wasteful, but may be worth 
investigating, as it would be simpler to implement. 
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7. The CVBEM was developed to solve 2D fluid flow problems fl^Q where 
components of the complex potential (the fluid potential $ or the streamline 
function ^E*) are known at different points around the contour, typically from 
physical measurements.^ A modification of the CBIEM can convert it to 
become a solver for Neumann (and thence mixed) boundary value problems. 
In the Neumann boundary value problem, C/^, the derivative of U across T, 
is known instead of U. Use the Cauchy-Riemann equations to observe that 
C/y = (the tangential component of V). The boundary information can 
be used to construct an approximation to V, by integration of Vt around F, 
using a suitable zero point (adding in a constant): 

v{j{t))^ f vAi{t'))dt'. 

Jo 

The same collocation process previously used to approximate V can in this 
case be used to approximate U. 

Beyond this, the technique is particularly applicable to free boundary prob- 
lems |8|, and may be able to be generalised to other elliptic (and possibly 
other second order) operators. 

8. The method would easily parallelise. The establishment of the linear system is 
computationally expensive, more so for high order interpolatory polynomials. 
This, as well as solution of the linear system, would efficiently (geometrically) 
parallelise. 



22 Other references to the CVBEM include [|l^, Fi, |l^ |^ . 

7 + zv is equiva 



23Warning: the notation used here W = U + j v^is equivalent to the fluid flow notation W = 
3> + i'i , so that U and V here have a different meaning from the fluid flow case, where they are 

commonly the components of the velocity q = {/i + Vj , and U = = , V = = . 

dx dy dy dx 
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A Listing of matlab ".m" files 



A.l cbiem.m 



function Vnnorm = cbiemCCCase, D, sigma, 0, alpha); 

■/.function Vnnorm = cbiemCCCase, D, sigma, 0, alpha); 
■/. 

■/, Perform the CBIEM on the Dirichlet problem. 

■/. 

■/. David De Wit March 30 1992 - January 14 1993 



■/. 0.0: Input and other parameters, togther with definitions. 



if "exist ( 'alpha' ) , 

if "exist('O'), 

if "exist (' Sigma' ) , 

if "exist('D'), 

if "exist ( 'CCase' ) , 



alpha = 2/3; end 
= 12; end 
Sigma = 0.32; end 
D = 7; end 
CCase = 4; end 



format short e; format compact; 

tvixixixixixixixixixixixixixixm^^^^^^ 
t 

'/, CCase: Index number to the contour being used: 



■/, 1: Unit circle, with 4 equally spaced artificial corners. 

■/, 2: Chandler's Teardrop, one corner, right angled, at the origin. 

■/. (This has been reversed to make it ACW. How different 

% from both GAC and DDW thesis.) 

■/. 3: Kress' ACW Teardrop, one corner, 2 pi/3 - angled. 

■/, 4: Kress' Reentrant contour, 3 pi/2 - angled. Reversed to 

■/, avoid the branch cut on the negative real axis, and 

■/. make it ACW in orientation. 

% 5: ACW Cardioid. Reentrant contour with 2 pi interior angle. 

% 6: ACW Heart. Reentrant contour with 2 pi and interior angles. 

% 7: ACW Controlled Cardioid, using trig parameterisation. 

■/, 8: ACW Controlled Cardioid, using polynomial parameterisation. 

■/. 9: Modified Boomerang, with a 5 degree external angle. 

•/. 

■/, D: Density of the geometric mesh, a positive integer. Choice 

■/, of D forces N, the size of the linear system being solved, 

■/. to be N = NC.(D+l)-2. 
■/. 

■/. sigma: Mesh parameter. Ratio of distances of consecutive mesh points 

■/, from the nearest comer. < sigma < 0.5. Try sigma = 0.25 

■/. as a starting guess. 
■/. 

■/, 0: Order of the interpolatory polynomial used to aproximate 

■/. V at the collocation points. Actually the (even) number of 

% nearest points used, thus 2 is linear, 4 is cubic. Must be 

% kept 2 <= <= NS+1; in practice keep <= 20. This limits 

% Q <= 8 for D = 2, 16 for D = 3, etc. 
•/. 

'/, alpha: Exponent of the true solution of test problem, W = z"{alpha}. 
•/. 

■/, NC: Numbers of corners/segments of the contour. 
7. 

'/, NS: Number of points on each side of the contour. N node points 

% create NS = (N+1) "2+1 mesh points on each segment . 

•/. 



mx/x/x/x/x/x/x/x/x/x/x/x/xm^ 
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1.1: Generate quadrature rule for a geometrically graded mesh, on 
one segment in the t domain using hprmesh, then insert this 
% into the grid with corners, using rmesh. 

numcorn =[4 111111111]; NC= numcorn(CCase) ; 



% graded mesh 
G = [0 Sigma. -(D:-l:l)]>; 
S = [2:D+2, D+l:-l:2] ' ; 
[t, w] = hprmeshCG, S, 0) ; 
N = length (tn); 



G = [G; 1-G(D+1:-1:1)] ; 

c = [0:NC] VNC; 

[tn, wn] = rmeshCc, t, w, 1) ; 



'/, Uniform mesh 

■/.%•/. N = (D+l)-2; tn = [1:N]'/N; wn = ones(size(tn))/N; 



"L 1.2: Compute the (complex) values of zn and zc (z at the node and 

% collocation points specified by tn and tc, respectively). 

% zc are found as the midpoints of zn in an arc-length sense, by 

■/o mapping the midpoints of tc to the contour. Also calculate 

% gdot, which is used to modify wn. 



tc = tn - diff([0; tn])/2; 

ptn = pi*tn; ptc = pi*tc; 

if (CCase == 1) 

zn = exp(2*i*ptn) ; zc = exp(2*i*ptc) ; 

gdot = 2*pi*i*exp(2*i*ptn) ; 
elseif (CCase == 2) '/, Chandler's Teardrop 

zn = 2*sin(ptn) - i*sin(2*ptn) ; 

zc = 2*sin(ptc) - i*sin(2*ptc) ; 

gdot = 2*pi*(cos(ptn) - i*cos(2*ptn)) ; 

gdot(N) = -2*pi*i; 
elseif (CCase == 3) '/, Kress 's Teardrop 

zn = sin(ptn) *2/sqrt (3) - i * sin(2*ptn) ; 

zc = sin(ptc) *2/sqrt (3) - i * sin(2*ptc) ; 

gdot = 2*pi*( cos(ptn)/sqrt(3) - i*cos(2*ptn) ); 

gdot(N) = -2*pi*i; 
elseif (CCase == 4) '/, Kress 's Boomerang 

a = 2/3; 

zn = - a * sin(3*ptn) - i * sin(2*ptn) ; 
zc = - a * sin(3*ptc) - i * sin(2*ptc) ; 
gdot = - 2*pi*( (3*a/2)*cos(3*ptn) + i*cos(2*ptn) ); 
gdot(N) = -2*pi*i; 
elseif (CCase == 5) '/, Plain Cardioid 

zn = (-1 + cos(2*ptn)) .*exp(i*2*ptn) ; 
zc = (-1 + cos(2*ptc)) .*exp(i*2*ptc) ; 

gdot = -2*pi*exp(i*2*ptn) .*(sin(2*ptn) + i*(l-cos(2*ptn))) ; 
gdot(N) = 0; 
elseif (CCase == 6) '/. Pointed Heart - Silly 

zn = - sin(3*ptn) - i*(sin(2*ptn) ) . "3; 
zc = - sin(3*ptc) - i*(sin(2*ptc)) ."3; 

gdot = - 3*pi*( cos(3*ptn) + i*2*cos(2*ptn) . *(sin(2*ptn) ) . "2 ); 
gdot (N/2) = 0; gdot(N) = 0; 

elseif (CCase == 7) 7. My Cardioid 

zn = - sin(3*ptn) - 5 * i* tn .* (1 - tn) .* sin(2*ptn) ; 
zc = - sin(3*ptc) - 5 * i* tc .* (1 - tc) .* sin(2*ptc) ; 
gdot = - 3 * pi * cos(3*ptn) - 5 * i * ... 

( (l-2*tn) .*sin(2*ptn) + 2*ptn.*(l-tn) .*cos(2*ptn) ); 
gdot(N) = 0; 
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elseif (CCase == 8) 7, My Stupid Polynomial Cardioid 
a = 7/24; 

zn = tn .* (tn - a) .* (tn - 1 + a) .* (tn - 1) + ... 

i * tn.'"2 .* (tn - 1/2) .* (tn - 1)."2; 
zc = tc .* (tc - a) .* (tc - 1 + a) .* (tc - 1) + ... 

i * tc."2 .* (tc - 1/2) .* (tc - 1)."2; 
gdot = (tn - a) .* (tn - 1 + a) .* (tn - 1) + . . . 

tn .* (tn-l+a) .* (tn-1) + tn .* (tn-a) .* (tn-1) + ... 

tn . * (tn-a) . * (tn-l+a) + . . . 

1 * ( 2 * tn .* (tn - 1/2) .* (tn - 1) . "2 + . . . 

2 * tn.*2 .* (tn - 1/2) .* (tn - 1) + ... 
tn.-2 .* (tn - 1) ."2 ) ; 

gdot(N) = 0; 

zn = 100*zn; zc = 100*zc; gdot = 100*gdot; 
elseif (CCase == 9) '/, 5 degree external angle. 

deg = 5; 

a = 2/( 3 * taii(deg*pi/360)) ; 
zn = - a * sin(3*ptn) - i * sin(2*ptn) ; 
zc = - a * sin(3*ptc) - i * sin(2*ptc) ; 
gdot = - 2*pi*( (3*a/2)*cos(3*ptn) + i*cos(2*ptn) ); 
gdot(N) = -2*pi*i; 
elseif (CCase == 10) '/, 20 degree external angle, 
deg = 20; 

a = 2/( 3 * taii(deg*pi/360)) ; 

zn = - a * sin(3*ptn) - i * sin(2*ptn) ; 

zc = - a * sin(3*ptc) - i * sin(2*ptc); 

gdot = - 2*pi*( (3*a/2)*cos(3*ptn) + i*cos(2*ptn) ); 

gdot(N) = -2*pi*i; 

end 

■/.rzn = real (zn) ; izn = Imag (zn) ; 

■/.plot (rzn , izn , ,rzn,izn, ' + ') ; grid; 

%plot(tn(l:10),izn(l:10) , ,tn(l:10) ,izn(l:10) ,'+'); grid; 
wn = wn.*gdot; 

■/, 2.1: Set up the complex order N matrices A and B. Establish F, a 
■/, matrix of indices to be used in calculating B, using Ft, a 

■/, submatrix of the pattern of indices for one edge. Also compute 

■/o L, the matrix of the coefficients of the interpolatory 

■/, polynomial, using the indices contained in F. 

oN = ones(N,l); j = 0/2; 

A = oN*wn.' ./ ( oN*zn.' - zc*oN'); NS = K/NC; 

Ft = zeros(NS, 0); Ft(l,:) = 0:0-1; 

Ft(NS, :) = NS-0+l:NS; 

Ft(2:j,:) = ones( j-1 , l)*Ft (1 , : ) ; 

Ft(NS-j+l:NS-l,:) = ones(j-l, l)*Ft (NS, : ) ; 

for k = j+l:NS-j, Ft(k,:) = Ft(k-1,:) + 1; end 

for k = 1:NC, F( (k-l)*NS+l :k*NS, : ) = Ft + (k-l)*MS; end 

F(l:j,l) = N*ones(j,l); 
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L = ones (size (F) ) ; ol = l:j; o2 = j+l:N; 

for k = 1:0, for v = 1:0 
if (v -= k) 

tk = tn(F(o2,k)); tv = tn(F(o2,v)); 

L(o2, k) = L(o2, k).* (tc(o2) - tv) ./ (tk - tv) ; 

tk = tn(F(l: j ,k)) ; tv = tn(F(l : j , v) ) ; 

if (k == 1), tk = 0; end; if (v==l) , tv = 0; end; 

L(ol, k) = L(ol, k) .* ( tc(ol) - tv ) ./ ( tk - tv ); 

end 

end , end 

IVXIX/XIXIX/XIXIXIX/W^ 

■/. 2.2: Compute the matrix B. J is a shift vector, used to create rows 

of B from rows of L. Multiply the temporary result through by 
% H, then write out C. 

H = suiii(A.'); tl = ones(j,l); 

t2 = [-tl; [0:NS-0-l]'; (NS-0)*tl]; J = zeros (size (H)) ; 

for k = 1:NC, J( (k-1) *NS+1 :k*NS) = t2 + (k-l)*NS; end 

B(l:j, 1:0-1) =L(l:j,2:0); B(l:j,N) =L(l:j,l); 

for k = j+l:N, B(k, J(k)+1 : J(k)+0) = L(k,:); end 

C = B.*(H'*ones(size(H))) - A; 

tvixixixixixixixixixixixixixixm^^^^^^ 

'I, 2.3: Set up and solve for Vn (computed approximation to V at the 
■/. node points), the linear system: 

■/. C Vn = A * Utn - sum(A) .*Uc 

Given C, establish the real, order N-1 matrix, Cstar, then 
■/, calculate Vn, using Vn(N) = 0. 

■/.sprintf (' Solving the linear system of size '/.g', N) 

Un = real(funccb(zn, alpha)); Uc = real(fmiccb(zc, alpha)); 

d = diag(H)*Uc - A*Un; d = -imag(d(2:N) - d(l:N-l)); 

C = real(C(2:N,l:N-l) - C(1:N-1,1:N-1)) ; 

Vn = zeros (size (Un) ) ; Vn(l:N-l) = C \ d; 

tixixixixixixixixixixixixixm^^ 

% 2.4: Error norms. Vne and Vce are the differences between the true 

and the computed values of V at the node and collocation points 
respectively, r is the discretisation error in the CBIE. p is 
the residual in the (above) computation for Vtn, when the true 

% soln at the node points is substituted into the equations. 

Vne = imag(funccb(zn, alpha)) - Vn; 
Vnnorm = sqrt (abs (wn) ' * (Vne . *2) ) ; 

■/. Compute the approximation at the 4 points of Kress: 
Wn = Un + i * Vn; 

zl = [ 0.1+0*i 0.2+0*i 0.3+0*i 0+0. 2*i ]; 
for j = 1:4 

WK(j) = sum(Wn .* wn ./ (zn - zl(j))) / sum(wn./(zn - zl(j))); 

end 

WKt = funccb(zl, alpha); Ek = [N abs(real(WK - WKt))] 
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A. 2 cint.m 



function [rho, C] = cint(sigma); 

7, function [rho, C] = cint (sigma) ; 
7. 

7, Contour integration with a geometric h-p grid. 
7. 

7, David De Wit July 14 1992 - July 17 1992 

%tt%ti:i:i:i:i:i:i:i:i:im^^^ 

if "exist ( 'QCase' ) , QCase = 2; end 
if "exist (' sigma' ) , sigma = 0.15; end 

DMin = 8; DMax = 15; 

format short e; format compact; 

QCase = 2; 

for D = DMin: DMax 

G = [0 sigma.-(D:-l:l)]'; G = [G; 1-G(D+1 : -1 : 1)] ; 

S = [2:D+2, D+l:-l:2]'; NC = 2; 

[t, w] = hprmeshCG, S, QCase, 0); c = [0:NC]'/NC; 

[tn, wn] = rmesh(c, t, w, 1); zn = exp(2*pi*i*tn) ; 

iN(D-DMin+l) = length(zn) ; 

gdot = 2*pi*i*zn; wn = wn.*gdot; 
ier(D-DMin+l) = 1-wn. '*( ( (zn-l)/i) . " (1/2) . /zn)/(2*pi*i*sqrt (i) ) ; 

end 

lier = loglO(abs(ier)) ' ; sqiN = sqrt(iN)'; 

plot(sqiN,lier, '-g' ,sqiN,lier, '+r') ; grid; 
title( 'Contour integration on a h-p geometric grid'); 

xlabel('sqrt(N) ') ; ylabeK 'loglO(Error) ' ) ; 
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A. 3 funccb.m 



function H = funccb(z, alpha); 

y, function W = funccb(z, alpha); 
% 

% True solution to the Dirichlet problem solved by cbiem. 

•/. 

•/, David De Wit April 13 1992 - September 27 1992 

%tt%ti:i:i:i:i:i:i:i:i:im^^^ 

W = z." (alpha) ; 



A. 4 funcci.m 



function W = f uncci(z) ; 

■/. function W = funcci(z); 
■/. 

Integrand of the problem solved by cint. 

■/. 

% David De Wit July 9 1992 - July 17 1992 



W = ((z-l)/i).-(l/2)./z; 



A. 5 funchp.m 

function f = funchp(x) 

■/. function f = funchp(x); 
■/. 

Function being integrated by hpmeth. 

■/. 

■/. David De Wit July 9 1992 - July 17 1992 



f = 1 - 3/2*sqrt(x); 
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A. 6 gettw.m 



function [QRt, QRw] = gettw(R, QCase) ; 

y, function [QRt, QRw] = gettw(R, QCase); 
% 

Get tables of nodes and weights for quadrature rules. User 
% inputs maximum number of points required, and the type required. 

The default type is Gauss — Lobatto (QCase = 2) , as it is of 
higher order than Newton — Cotes. 

■/. 

•/. David De Wit July 13 1992 - July 17 1992 

tttvi:i:i:i:i:i:i:i:i:i:i:im^^^ 

if "exist('R'), R = 20; end 
if "exist ( 'QCase' ) , QCase = 2; end 

tttvi:i:i:i:i:i:i:i:i:i:i:im^^^ 



if (QCase == 1) 

if (R > 10) 



R = 10; 

sprintfCR too large 



for Newton — Cotes. Now R = 10. \n'); 



end 

QRw = [ 



1 


1 


1 


7 


19 


41 


751 


989 


2857 


16067 


1 


4 


3 


32 


75 


216 


3577 


5888 


15741 


106300 





1 


3 


12 


50 


27 


1323 


-928 


1080 


-48525 








1 


32 


50 


272 


2989 


10496 


19344 


272400 











7 


75 


27 


2989 


-4540 


5778 


-260550 














19 


216 


1323 


10496 


5778 


427368 

















41 


3577 


-928 


19344 


-260550 




















751 


5888 


1080 


272400 























989 


15741 


-48525 


























2857 


106300 





























16067 



] 

QRw = QRw(l:R+l,l:R) ; 

for j = 1:R 

QRw(l:j+l,j) = QRw(l 
QRt(l:j+l,j) = [0:j]^ 

end 

elseif (QCase == 2) 

QRt = zeros (R+1,R) ; 
for j = 2:R+1 

[QRt(l: j ,j-l) , QRw(l 



QRt = zeros (QRw) ; 



j+l,j)/sum(QRw(: , j)) ; 

/j; 



QRw = QRt; 
j,j-l)] = lobatto (j ,0,1) ; 



end 



end 
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A. 7 hpmeth.m 



function [lEip] = hpmeth(DMax, GMax, p, sigma, QCase) ; 

y, function [lEip] = hpmeth(DMax, GMax, p, sigma, QCase); 
% 

% Experiment with h-p integration methods. 

■/. 

•/, David De Wit July 9 1992 - July 17 1992 

%tt%ti:i:i:i:i:i:i:i:i:im^^^ 

if "exist (' QCase' ) , QCase = 2; end 

if "exist (' sigma' ) , sigma = 0.15; end 

if "exist Cp'), p = 6; end 

if "exist (' GMax' ) , GMax = 6; end 

if "exist ( 'DMax' ) , DMax = 19; end 

tttvi:i:i:i:i:i:i:i:i:i:i:im^ 

'I, Geometrically-graded h-p method. 

Ehp = zeros (DMax, 1) ; vN = Ehp; 

for D = l:DMax 

G = [0 Sigma. -(D:-1:0)] ' ; S = [2:D+2] ' ; 

[t, w] = hprmeshCG, S, QCase, 0); 

Ehp(D) = funchp(t) '*w; vNhp(D) = length (t ) ; 

end 

IvNhp = loglO(vNhp); lEhp = loglO(Ehp); 

ituT/x/x/x/x/x/x/x/x/xm 

'I, Various linear h and p methods. Variable g, constant p. 

Eip = zeros (DMax, GMax) ; 
for D = l:DMax 
N = 2*D; 
for g = l:GMax 

G = ([0:N]'/N).-g; 

[t, w] = hprmeshCG, p, QCase, 0); 

vN(D) = length(t); Eip(D,g) = f michp(t) ' *w; 

end 

end 

IvN = loglO(vN); lEip = loglO(Eip); 

sprintf (' Order '/.f, slopes of the blue lines are approximately', p) 
(lEip(DMax-l, :) - lEip(DMax, :)) ./(IvN(DMax-l)-lvN(DMax)) 

plot(lvNhp,lEhp, '+r' ,lvNhp,lEhp, '-g' ,lvN,lEip, '+r ' , IvN, lEip, ' -b' ) ; 
xlabelCloglO (points) ') ; ylabel (' loglO (error) ') ; 

text(0.9,0.7, 'g = l','sc'); text (0 . 9 ,0 . 57 , ' g = 2','sc'); 

text(0.9,0.45, 'g = 3','sc'); text(0.9,0.35, 'g = 4','sc'); 

text(0.9,0.27, 'g = 5','sc'); text(0.9,0.2, 'g = 6','sc'); 

text(0.7,0. 15, 'h-p method' , 'sc') ; grid; 

ttTi:i:i:i:i:i:i:i:i:i:i:i:m^^^^^^ 
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■/. Compare h methods for quadrature rules of various p (# points) . 



OMax = 12; 

for p = 2: 2: OMax 

g = p; j = p/2 + 1; 

for D = l:DMax 

N = 2*D; G = ([0:N] VN) ."g; 

[t, w] = hprmeshCG, j, QCase, 0); 

Ehhp(D,j-l) = funchp(t)'*w; hliN(D) = length (t ) ; 

end 

end 



IhhN = loglO(hhN); lEhhp = loglO(Ehhp); 

sprintf (' Slopes of the blue lines are approximately') 
(lEhhp(DMax-l, :) - lEhhp(DMax, :)) ./(IhhN(DMax-l)-lhhN(DMax)) 



plot (IvNhp , lEhp , '+r' ,lvNhp,lEhp, '-| 
xlabeK'loglO (points) ') ; 
text(0.9,0.76, 'p = l','sc'); 
text (0.9,0. 5, 'p = 5','sc'); 
text(0.9,0.32, 'p = 9','sc'); 
text(0.7,0. 15, 'h-p method' , 'sc') ; 



' ,lhhN,lEhhp, ' +r ' , IhhN , lEhhp , '-b') ; 
ylabel ( ' loglO (error) ' ) ; 
text(0.9,0.62, 'p = 3','sc'); 
text(0.9,0.4, 'p = 7','sc'); 
text(0.85,0.2, 'p = 11', 'sc'); 
grid; 
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A. 8 hprmesh.m 



function [t, w] = hprmesh(G, S, IClosed) 
y, function [t , w] = hprmesh(G, S, IClosed) 



% 

Create a new quadrature rule, based on a mesh G, where between 

points G(i) and G(i+1) is a (closed) quadrature rule of 

Gauss — Lobatto type on S(i) points, including the 2 end points 

■/, G(i) and G(i+1); for i = 1 : length(G) -1 . If IClosed is 1, then 

the contour is closed, and the ends are tied together. This 

'/o function is a generalisation of rmesh. 

•/. 

■/. David De Wit July 13 1992 - December 2 1992 



if "exist('G'), sigma = 0.1; G = [0 sigma. " (3: -1 : 1) 1]'; end 

if "exist CS'), S = [2:5]'; end 

if "exist (' IClosed' ) , IClosed = 1; end 



IG = length(G); IS = length(S); 

if ((IG "= lS+1) k (IS "= D) 

sprintf Chprmesh: Danger KG) = '/,f, 1(S) = 'IS' , IG, IS) 

end 

Set up S for rules with a constant integration rule, 
if (IS == 1), S = ones(lG-l,l)*S; end 



dG = diff(G); S = S - 1; 

N = max(S); R = length(dG); 

7. Obtain the nodes and weights in a table 

QRt = zeros(N+l,N) ; QRw = QRt; 

for j = 2:N+1 

[qRt(l:j,j-l), QRw(l:j,j-l)] = lobattoCj ,0, 1) ; 

end 



% Play with the table 

QRw(N+l,:) = diag(qRw(2:N+l,l:N)) ' ; QRtCN+l,:) =ones(l,H); 

for 1 = 2:N, for j = l:i-l, qRt(i,j) = NaN; QRw(i,j) = NaN; end, end 

tt = QRt(:,S); tw = QRw(:,S); 

j = ones(N+i,i); tw = tw.*(j*dG'); 
tt = tt.*(j*dG') + j*G(l:R)'; 
tw(l,2:R) = tw(l,2:R) + tw(N+l , 1 :R-1) ; 

twl = tw(l:M,:); ttl = tt(l:N,:); 

t = [ttl(:); tt(N+l,R)]; t(isnan(t)) = [] ; 

w= [twKO; tw(N+l,R)]; w(isnan(w)) = []; 



if (IClosed == 1) 

N = length (t); 

w = [w(2:N-l); w(l)+w(N)]; t = t(2:N); 

end 
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A. 9 lobatto.m 



function [x, w] = lobatto(n, a, b) 

y, function [x, w] = lobatto(n, a, b) 
% 

Return the weights w and points x of the n-point Gauss — Lobatto 
% quadrature rule on the interval [a, b] . 

■/. See G. H. Golub, SI AM Review 1973 p 318. 

■/. 

•/, Graeme Chandler July 1992 

n = round (n) ; 
if (n == 2) 

X = [a; b]; w = [1; l]*(b-a)/2; 

elseif (n == 3) 

X = [a; a+(b-a)/2; b] ; w = [1; 4; l]*(b-a)/6; 

(n >= 4) 

nn = n-1; m = l:2:2*mi-l; 

m = (l:nn-l) ./ sqrt (m ( 1 : nn- 1 ) .* m(2:nn)); 
J = (diag(m,-l)+diag(ni,l)) ; 

1 = eye(nn); en = (l:nn)' == nn; 

gam = (J + I)\en; mu = (J - I)\en; 

sol = [1 -gain(nn); 1 -mu(mi)] \ [-1 ; 1]; 

alpha = sol(l); beta = sqrt (sol(2) ) ; 

[ww.xx] = eig([J beta*en; beta*en' alpha]); 
[xx, i] = sort(diag(xx)) ; 
w = ww(l,i) ' .~2 * (b-a) ; 
x= [a; (a+b)/2+(b-a)*xx(2:nn)/2; b] ; 
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A. 10 rmesh.m 



function [t, w] = rmesh(G, tl, wl, IClosed) ; 

y, function [t , w] = rmesh(G, tl, wl , IClosed); 
% 

Create a new quadrature rule, based on a mesh G, where a 
(closed) quadrature rule (tl, wl) is inserted over each 
interval of G. If IClosed is 1, then the contour is closed, 
and the ends are tied together. This function is generalised 
into hprmesh. Originally conceived by Graeme Chandler. 

■/. 

% David De Wit July 13 1992 - September 6 1992 

•/. 

if (length(G) < 2), return; end 

if (length (G) "= 2) 

h = diff(G); tw = wl*h' ; 

n = length (h) ; m = length (wl ) ; 

tw(l,2:n) = tw(l,2:n) + tw(m,l:n-l); twl = tw(l :m-l , : ) ; 

w = [twl(:); tw(m,n)] ; 

tt = tl(l:m-l)*h' + ones(m-l,l)*G(l:n) ' ; 
t = [tt(:); G(n+1)]; 

else 

t = tl; w = wl; 

end 

if (IClosed == 1) 

N = length(t) ; 

w = [w(2:N-l); w(l)+w(N)]; t = t(2:N); 

end 
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A. 11 testcb.m 

function N = testcb(C, sigma, Dmin, Dmax, Omin, Omax, alpha) 

y, function N = testcb(C, sigma, Dmin, Dmax, Omin, Omax, alpha) 
% 

% Run cbiem for various parameters, and tabulate results. 

■/. 

•/, David De Wit May 12 1992 - December 21 1992 

%tt%ti:i:i:i:i:i:i:i:i:im^^^ 

if "exist ( 'alpha' ) , alpha = 1/2; end 
if "exist ( 'Omax' ) , Omax = 16; end 
if "exist (' Omin' ) , Omin = 8; end 
if "exist ( 'Dmax' ) , Dmax = 19; end 
if "exist('Dmin') , Dmin = 16; end 
if "exist (' sigma' ) , sigma = 0.32; end 
if "exist CO, C = 7; end 

format short e; format compact 



for D = Dmin: Dmax 

for = Omin: 2: Dmax 

P = cbiem (C, D, sigma, 0, alpha); 
if ((P -= Inf) & (P "= NaN)) 

N(D-Dmin+1, (0-Omin) /2+1) = P; 

else 

N (D-Dmin+1 : Dmax-Dmin+1 , : ) = ... 

Inf *ones(Dmax-D+l, (0max-0)/2 + 1) ; 

return 

end 

end 

if (min(N(D-Dmin+l, :)) > 1) 

H(D-Dmin+2:Dmax-Dmin+l, :) = ... 

Inf *ones(Dmax-D, (0max-0min)/2+l) ; 

return 

end 

end 
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